Theory of rigid-plane phonon modes in layered crystals 
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The lattice dynamics of low-frequency rigid-plane modes in metallic (graphene 
multilayers, GML) and in insulating (hexagonal boron-nitride multilayers, BNML) 
layered crystals is investigated. The frequencies of shearing and compression (stretch- 
ing) modes depend on the layer number N and are presented in the form of fan 
diagrams. The results for GML and BNML are very similar. In both cases only 
the interactions (van der Waals and Coulomb) between nearest-neighbor planes are 
effective, while the interactions between more distant planes are screened. A com- 
parison with recent Raman scattering results on low-frequency shear modes in GML 
[Tan et al, larXiv: 1106. 1146V 1 (2011)] is made. Relations with the low-lying rigid- 
plane phonon dispersions in the bulk materials are established. Master curves which 
connect the fan diagram frequencies for any given ?\f are derived. Static and dynamic 
thermal correlation functions for rigid-layer shear and compression modes are calcu- 
lated. The results might be of use for the interpretation of friction force experiments 
on multilayer crystals. 



INTRODUCTION 



The experimental discovery of graphene and other free-standing two-dimensional (2D) 
crystals [l|, |2| has opened the path for the synthesis of a whole class of layered materials with 



2| (3D h-BN 
6|]). While graphene is a purely covalent 



novel physical properties and with a great potential for technological applications jsl- The 
most prominent member, graphene — a monoatomic layer of crystalline C with hexagonal 
structure — is a metallic conductor. In addition to an unusual electronic spectrum, this 
material shows extraordinary mechanical strength [4] and thermal properties 5|. 

On the other hand 2D hexagonal boron-nitride (h-BN) is an insulator jj, 
has a direct band gap in the ultraviolet region 
crystal, 2D h-BN, built from III-V elements, has partially covalent and ionic bonds. The 
ionic character is a consequence of the charge transfer of ~ 0.6 electrons from B to N 
Since the crystal structure of 2D h-BN is non-centrosymmetric (point group symmetry D^h), 
the two sublattices (B"*" and N~) exhibit an electromechanical coupling. Hence 20 h-BN is 
the structurally most simple crystal which, according to theoretical predictions [8], should 
be piezoelectric. 

Nanoscale thin sheets of graphene, 2D h-BN and related layered materials Q| are of great 
importance for applications as electronic devices and nanoelectromechanical systems. The 
synthesis and characterization of multilayers and the study of their physical and chemical 
properties is a challenge of current solid-state physics and materials science. In particular, 
the change of properties with the number of layers and the evolution of the layer system 
to the corresponding bulk material are of foremost importance. Most remarkable is for 
instance the change in electronic structure from graphene, a zero-gap semiconductor, to 
graphite, a semimetal with band overlap [l^. These theoretical results are directly related 



to the interpretation of electronic transport experiments The change in the 



electronic bands is reflected in the double resonance Raman spectrum that clearly evolves 
with the number of layers 13]. Beside the electronic structure, the elastic properties depend 
on the number of layers. Atomic force microscopy experiments (AFM) on various thin- 
sheet materials demonstrate that the nanoscale friction decreases with increasing number of 
layers It has been suggested that the trend arises from the thinner sheets' increased 
susceptibility to out-of-plane elastic deformations. 

Since physical properties vary with the number of layers, it is important to study sep- 
arately the lattice dynamics of modes where the atomic planes move as rigid units. Since 
this motion is governed by the weak interlayer forces, the corresponding frequencies are low 
(< 150 cm~^) in comparison with the high-frequency optical modes (< 1600 cm~^) which 
are due to covalent intralayer forces. In graphite these low-frequency modes have been 



discovered half a century ago by inelastic neutron scattering experiments [l^. One distin- 
guishes modes where the planes move parallel to the hexagonal axis (we will call these modes 
compression modes), and modes where the planes move perpendicular to this axis (shear 
modes). Later on the complete phonon dispersions associated with rigid-plane motion have 



been measured by neutron scattering [15!]. The rigid-plane (-layer) shear mode is optical 
active and has been measured by Raman scattering in graphite [ig] and in 3D h-BN (l7| . 
Due to the weakness of the interlayer forces, the rigid-layer shear frequency in graphite [18 1 



and in h-BN [19| increases strongly with applied pressure. 

The measurement of rigid-layer modes in few-layer systems has been an outstanding 
problem. Neutron scattering is not an adequate technique since the samples are too small. 
Most recently, the interlayer shear modes in few-layer graphene systems have been uncovered 



by Raman spectroscopy 20|]. The increase of the resonance frequency with increasing layer 
number provides a unique signature for few-layer graphene systems and for multilayers in 
general. 

In the present paper we report on theoretical studies of low-frequency rigid-layer shear 
modes and compression modes in graphene- and boron-nitride multilayers. While the high- 
frequency optical mode spectra of graphene- and boron-nitride multilayers are very different 
due to the efficiency of Coulomb forces in the latter |21], the low-frequency optical mode 
spectra in both systems turn out to be very similar. 

The content of the paper is as follows. First (Sect. II) we present the main theoretical 
concepts which are used to treat by analytical means the lattice dynamics of multilayer 
systems. Within a same formalism we consider graphene multilayers and h-BN multilayers. 
The phonon dispersion relations of the corresponding rigid-layer motions (compression and 
shear modes) are calculated in Sect. III. The dependence of the frequency spectra on the 
number of layers K is presented in the form of fan diagrams. Next (Sect. IV) we compare 
the theoretical results with experiment. Then we derive relations between the phonon fre- 
quencies of the rigid layer systems and the dispersions of the corresponding bulk materials. 
We derive master curves which allow to connect the fan diagram frequencies for any given 
K. In Sect. V we calculate static and dynamic thermal displacement correlation functions. 
The temperature dependence of displacement correlations of surface layers is calculated, the 
dependence on the layer number is investigated. Concluding remarks (Sect. VI) close the 
paper. 



FIG. 1: Brillouin zone of the 3D hexagonal primitive lattice F/^; the shaded hexagon containing 
the F, M and K points is the Brillouin zone of a 2D hexagonal crystal. 



II. LATTICE DYNAMICS 



In previous work we have studied by analytical methods the phonon dispersion relations 
for graphene multilayers (GML) 22| and h-BN multilayers (BNML) 2lJ. Here we briefly 
recall the main concepts. Both graphene and 2D h-BN have the same symmetry D^h with 
two atoms per unit cell. Since each atom has three degrees of freedom, the dynamical matrix 
for the planar problem has dimension 6x6. Here gl is the wave vector in the 2D 
Brillouin zone (Fig. [1]). The 3D parent crystals, graphite and bulk h-BN have the same 
space group symmetry, PG^/mmc (D^f^). Since in both cases there are 4 atoms per unit 



cell, the corresponding dynamica. 
electron diffraction experiments 
same as in graphite (. . .ABAB. 



matrices are of dimension 12 x 12. In the case of GML, 



12 



have shown that the stacking of atomic planes is the 



23| . In 3D h-BN, each B atom is on top of a N atom 



in the adjacent plane and vice versa, with . . .AA'AA' . . . stacking 2J|. We assume that the 
same holds for BNML. 

We will use a unified description of the dynamical matrix for GML and BNML; the 
differences in structure and in interatomic (ionic) forces will be taken into account in the 
numerical evaluation of the secular equation. We consider a slab of a finite number of 'N 
layers, the layers are labelled by an index / G {0, 1, . . . , Jsf — 1}. The distance between 
nearest-neighbor planes which are perpendicular to the crystallographic c axis is c/2. Since 
the slab is indefinitely extended only in two dimensions, we regard it as a 2D crystal which 
consists of prismatic unit cells 25| with basis area 73/2 and height Kc/2. We recall 
that a = 1^1 1 02! is the length of the lattice translation vectors of the 2D hexagonal 



basis crystal 



23| . Each unit cell contains 3\f pairs of atoms (C,C) or (B,N) in the case of 
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GML or BNML, respectively. Since each atom has three vibrational degrees of freedom i 
(j) G {x, y, z}, the !N-layer slab has 6IN" vibrational modes. 

In order to calculate the phonon dispersion relations, we construct the 63\fx 6K dynamical 
matrix Aj^(q±). In terms of 6 x 6 submatrices D(/,/'|gl) with elements D^j^' (l,l'\qj_), where 
K (k') takes two values k {k') G {C, C} or {B,N} which corresponds to {C, C} or {B,N}, 
the dynamical matrix reads 



/ 2)(0,0|g-l) 
25(1, 0|gl) 
2)(2,0|g-l) 



2)(0,l|gl) 
25(1, l|gl) 
2)(2,l|gl) 



2)(0,2|g-l) 
25(1, 2|gl) 

2)(2,2|g-l) 



l,0|gl) T'(K-l,2|gl) ... - 1, X - l|gl) / 



T>(0,X-l|gl) \ 
2)(2,X-l|gl) 



(1) 



Here we take into account the interaction within a same layer (/ = /') and interactions 
between layers separated by a distance (/ — l')c/2^ I ^ V . The "same-plane" matrices 
2^(^5^|^l) ^-re given by 



D(/, = D{1, /Igl) + K{1, l\u = 0), 



(2) 



where D{1, l\q±) is the dynamical matrix of the l-th monolayer while K{1, l\q± = 0) accounts 
for the self-interaction due to interplane couplings. Assuming that the in-plane interactions 
are the same for all planes, one has in terms of elements 



D^fiUlq^) = Ft/{l,l\q^) + C^f{l,l\qi. 



(3) 



The matrices F and C stand for the in-plane covalent and Coulomb interactions, respectively. 
In the case of GML, only covalent interactions are taken into account, F 7^ 0, C = 0; in 
the case of BNML, both F and C are taken into account. The interplane coupling matrices 
D(/,/'|gl), / 7^ /' in Eq. ([1]) are due to van der Waals and Coulomb interactions J and C, 
respectively: 



Vtfil,l'\q±) = J^''{q^)5i',i±i + Qr'il^^'W±)- 



(4) 



Here again Coulomb interactions are relevant for BNML with 1<|/ — /'|<3\f— 1, while 
for both BNML and GML van der Waals forces act between nearest neighbor planes only. 
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Hence in the case of GML, only nearest-neighbor ofF-diagonal elements D{l,l ± l\q±) are 
non-zero in Eq. (P (see Eq. (27) of Ref. ^ ). 

In calculating the submatrices F{l,l\q±) we take into account intra-plane covalent inter- 
actions by means of a force-constant model originally derived from in-plane inelastic X-ray 
scattedug expedmeuts on single crystals of gxaphite Q. In case of BNML. intra- and 
extra-plane Coulomb interaction matrices are calculated by means of Ewald's method [27J. 
We have solved numerically the secular determinant of order 6!N, 



0, 



(5) 



and obtained the phonon dispersion relations for GML j22| and BNML 2l|. In terms of 
eigenvalues [u}\{q±)Y ^''^^ orthogonal eigenvectors ^x{q±) we have 



(6) 



W kk' ij 



Here A labels the 6IN" eigenmodes. We discern 3 acoustic modes such that uj\{q±_ = 0) = 
for A = 1, 2, 3 and 6K — 3 optical modes with UJ\{q±) ^ for all values of q±. Among the 
latter we distinguish near = 3?\f atomic vibrational modes (in-plane and out-of-plane 



displacements) and 3(?sf — 1) rigic. 
As has been shown previously 



ane modes. 

2l|, marked differences between GML and BNML appear 
in the highest optical branches with frequencies ~ 1300 - 1500 cm~^. These modes are due 
to intra-plane shear displacements where the Coulomb forces in BNML are efficient. On the 
other hand for a given 3\f, the low-frequency (< 200 cm~^) optical phonon dispersions in 



GML and BNML are very similar 2l|, |22|. These modes are due to rigid-plane compression 



and shear displacements. Due to the overall charge neutrality of the BN atomic planes, 
the Coulomb forces between nearest-neighbor planes are screened and hence there is no 
qualitative difference between GML and BNML rigid-plane modes. In the following section 
we will discuss the evolution of the rigid-layer modes with the number of layers. 



III. RIGID-LAYER MODES 



We first recall the situation in the 3D parent materials graphite and h-BN. Since the 
point group symmetry is -De/i, the decomposition into irreducible representations of the 
optical displacements at the center of the Brillouin zone (Fig. [1]) reads j28|, l29| F = + 
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2B2g+Eiu + 2E2g. While six of these modes are high-frequency inter-plane vibrational modes 
(800 - 1600 cm~^), one of the doubly degenerate modes and one mode refer to low- 
frequency rigid-plane motions. We denote these modes by £'2^1 and B2g^. The E2g^ mode 
corresponds to the rigid-plane shear displacements perpendicular to the crystallographic 
c axis. This mode has been measured in graphite by Raman reflectivity [29] at 42 ± 1 
cm~^ and is called "rigid-layer shear" mode. It can be identified with the zero wave vector 
transverse optical mode (TOj^near 1.35 THz measured first by neutron scattering in high- 
quality pyrolytic graphite [151]. The i?23i mode corresponds to a rigid-layer compression 
mode along the c axis. This mode is optically inactive, however it appears near 3.9 THz 
130 cm~^) in neutron scattering 15|]. For a more complete discussion of the early work, 
we refer to Ref. [16|]. A discussion of the zone-center optical modes in h-BN was originally 



given in Ref. 30|], the i?2gi rigid-layer shear mode was observed by Raman scattering 17| at 
51.8 cm~^. Recently the phonon dispersions of h-BN have been measured by inelastic x-ray 
scattering and analyzed by ab initio calculations 311] • At the F point, the i?2gi rigid-layer 
shear mode has energy 6.5 meV (52 cm~^), the compression rigid-layer mode (called there 
Big) has energy 15 meV (121 cm~^). 

In the !N-layer system the 3(3\f — 1) rigid-plane optical modes at the center of the 2D 
Brillouin zone decompose into !N — 1 compression modes with frequencies u)x^{q± = 0), 
Ac = 1, 2, 3, . . . , 3\f — 1 and eigenvectors ^(Ac, q± = 0) and into 3\f — 1 doubly degenerate 
shear modes with frequencies ux^{q± = 0) and eigenvectors ^(aI^"*, q± = 0) or ^{X^s'\ q± = 0), 
where A^^'' (A^^^) = 1, 2, 3, . . . , ?sf — 1. The degeneracy of the shear modes is a consequence of 
hexagonal symmetry at q± = 0. The compression and shear modes correspond respectively 
to the low-frequency modes -8231 and E2gi of the bulk materials. Since the Ac modes have 
only nonzero displacement components along the z-direction (c-axis) while the A^ modes 
have only nonzero x- and y-components, one has the simplified orthonormality conditions 

E^^^(^-0)^^^(^-0) = 5a.a^, (7) 

l,K 

5^5^ep)(A„0)ef''')(A;,0) = 5,,,,, (8) 

i l,K 

where i G {x,y}. The eigenvector components fulfill the relations 
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where m^^ is the mass of particle k, Aq, = Ac (A^) for i = z {x,y). Hence particles within a 
same plane / experience equal displacements. In addition one has 

J2^f'''\Xa,0)=0, (10) 

I 

where a = s for i G {x, y} and a = c for i = z, which means that the center of mass of the 
N-layer system stays at rest. 

As an example we first consider GML. For the case 'N = 2 (bilayer), we obtain the 
rigid-plane shear mode eigenvectors with 12 components 

6^=2(x) = ^(1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0), (11) 
ei=2(l/) = ^(0, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0), (12) 

with degenerate eigenfrequency 30.4 cm~^ [see Fig. [2t^a)]. Here the components 1 — 3 and 
4 — 6 refer to the Cartesian displacements of the first (k = 1) and second {k = 2) particle in 
plane / = 1, respectively, while the components 7 — 9 and 10 — 12 refer to the displacements 
of the two particles in the plane 1 = 2. The rigid-plane compression mode eigenvector for 
the bilayer reads 

6,=2(^) = ^(0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1), (13) 

with eigenfrequency 90.1 cm~^ [see Fig. |2]^b)]. 

For the case 3\f = 3 there are two doubly degenerate shear modes with frequencies 21.5 
cm~^ and 37.2 cm~^, shown respectively in Figs. M^c) and (e), and two compression modes 
with frequencies 63.7 cm^^ and 110.4 cm~^, Figs. [2](d) and (f). We notice that for one 
degenerate shear mode and one compression mode [Figs. [21(c), (d)] the center plane undergoes 
no displacement. This feature is characteristic for all multilayers with ?sf uneven. On the 
other hand the displacements of the two outer layers are opposite to the displacements of 
the inner (central) layer [Figs. [2](e) and (f)], in agreement with the general requirement 
that the center of mass of the !N"-layer system stays at rest. We have plotted the calculated 
frequencies of the rigid-plane shear modes and compression modes for GML as a function of 
the layer number K in the form of fan diagrams in Fig.[3]^a). The lower set (red, filled circles), 
centered around the bilayer shear frequency 30.4 cm~^, corresponds to the rigid-layer shear 
modes {co'As(^)}; the upper set (blue, open circles), centered around the bilayer compression 
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a) ^ b) 1^ 




FIG. 2: Rigid-plane shear and compression displacements, }sf = 2, (a) and (b), respectively; X = 3, 
(c), (e) and (d), (f), respectively. 

frequency 90.1 cm^^ corresponds to the rigid-layer compression modes {(X'a^(^)}. In the fan 
diagram of shear motions for GML the sequence of frequencies at 30.4 cm~^ which occurs 
for even values of !K refers to the situation where the ^ upper planes of the system move 
in unison in a direction perpendicular to the hexagonal axis (c), while the ^ lower planes 
move in unison in opposite direction. The same holds for the compression motions where 
the frequencies at 90.1 cm~^ in the fan diagram correspond to two sets of unison motions 
with opposite direction along c. 

Proceeding along the same lines, we have calculated the rigid-layer eigenfrequencies and 
displacement vectors for the case of BNML. Shown in Fig. [3t^b) are again the eigenfrequency 
fan diagrams as a function of 'N. The lower set (red, filled circles) corresponds to the rigid- 
layer shear modes and is centered around the BN bilayer shear frequency 38.6 cm~^, the 
upper set (blue, open circles), centered around the BN bilayer compression frequency 86.3 
cm^^ corresponds to the rigid-layer compression modes. 

In the next section we discuss these results and compare with recent experiments and ab 
initio calculations on GML 12011. 
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g 



g 




FIG. 3: Rigid-layer frequencies a;(!N) at the T point for (a) GML and (b) BNML as a function of 
number of layers ?\f. Lines connecting the frequency points are master curves given by Eqs. (jlSD 
and (jl9p . For both GML and BNML, the two sets of points (red, filled circles and blue, open 
circles) are referred to as lower and upper fan diagrams, respectively. 



IV. DISCUSSION OF RESULTS 



A. Comparison with experiment 

Given the similarity of the fan diagrams (Figs. [3t^a) and (b) for GML and BNML, re- 
spectively) we first discuss some general features. We recall that for the GML systems we 
have only taken into account van der Waals forces between nearest-neighbor planes. On 
the other hand we have treated the BNML systems as ionic crystals. In addition to van 
der Waals forces between nearest-neighbor planes we have summed the Coulomb interac- 
tions over all planes. Obviously the overall charge neutrality and the layer rigidity, i.e. 
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the absence of relative motion between the B"*" and N~ sublattices, leads to a screening 
of the Coulomb interactions between next-nearest-neighbor and more distant h-BN planes. 
Hence the low-frequency fan diagrams for GML and BNML are very similar. We recall that 
the screening effect disappears for non-rigid-layer displacements, as is seen from the high- 



frequency (> 1300 cm ^) optical mode dispersions in BNM 



'221. 



21| which are qualitatively 



very different from the corresponding dispersions in GML 

The similarity of the low-frequency results of GML and BNML is also reflected by the 
low-frequency spectra of the 3D parent materials. We recall that transverse acoustic (TA) 
and optical (TO) as well as longitudinal acoustic (LA) and optical (LO) phonon dispersions 
along q = (0, 0, Qz) were first measured by neutron scattering ISj in pyrolytic graphite and 
recently by inelastic X-ray scattering in single crystals of graphite 261] and in h-BN 3l| . In 
igs. 111(a) and (b) we ha ye p lotted the low-frequency phonon dispersion relations for graphite 



22l | and for bulk h-BN 2l| along the line A-F in the 3D Brillouin zone (see Fig. [T 



. Similar 



33| and 



results have been obtained earlier by first-principles calculations for graphite [32|, 
h-BN 

We first discuss the low-frequency dispersions of graphite in relation with the GML low- 
frequency shear modes. The two lowest branches (TA and TO) in Fig. 111(a) (red, full lines) 
refer to the acoustic and optical rigid-plane shear motion, respectively. The optical branch 
TO evolves from the frequency value 30.4 cm~^ at the A point of the Brillouin zone to 
43.0 cm~^ at the F point. We notice that the value at the A point [qa = (0, 0, ^)] agrees 
with the rigid-plane shear frequency 30.4 cm~^ of the bilayer at gl = and with the 
center points for even !N in the lower fan diagram of Fig. [31(a) (red, filled circles). We 
attribute this agreement to the fact that in our calculation on graphite and multilayer 
systems we have restricted ourselves to van der Waals interactions between nearest-neighbor 
planes only I22I]. However, since this is a well-justified approximation, this correspondence 
should also be valid experimentally (inelastic neutron scattering on graphite [l^ gives jy ^ 1 
THz at the A point, most recent low- frequency Raman experiments 20| measure a shear 
frequency of 31 cm^^ for bilayer graphene). The value 43.0 cm~^ of the £'291 niode at the 
F point of the 3D system corresponds to the limit value for large 3\f of the sequence of 
highest frequencies {u;^^([N")} = {30.4 cm~^; 37.2 cm~^; 39.7 cm~^; . . .; 42.4 cm~^; . . .} for 
K = {2; 3; 4; . . . ; 10; . . .}, respectively. 

In case of h-BN [Fig. IH^b)] the rigid-layer shear mode TO evolves from 38.6 cm~^ at 
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FIG. 4: Low-frequency phonon branches along A-F {q = {0,0, q^) : (0,0, ^) — > 0) for (a) graphite 
and (b) 3D h-BN. Cuts at well-defined values of qz (marked by vertical lines and labels at the top 
horizontal axis) yield the F point frequencies of the multilayers (see text). 

A to 54.5 cm"^ at F (symmetry E2g^). Notice that the value at A agrees again with the 
shear mode frequency of the bilayer. Here the Coulomb interaction is only effective between 
nearest-neighbor layers and screened between more distant layers in the bulk system and in 
[N"- layer systems. The value at F agrees with the corresponding limit frequencies {a;^^(!N')} 
of BNML for ?^ — ^ oo. 

The sequence of frequencies {u^ (Hi) been measured recently in GML up to K = 



11 by polarized Raman techniques 



20|. These experiments demonstrate that the shift 



of the resonance (called C peak |20 [) with !N is truly representative of the GML system. 



Furthermore the authors of Ref. 



20[ have studied the eigenfrequencies and eigenvectors 



of the 3\f-layer shear modes by using a simple linear-chain model and by performing ab 
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initio calculations up to [N" = 5. Comparing our own calculated values of the shear- mode 



eigenfrequencies with those of Ref. 



20|, we see close agreement. 



It is useful to trace back this agreement on the level of interlayer van der Waals (vdW) 
interactions. From experiment the authors of Ref. 2o| derive that the interlayer force constant 
per unit area a has the value a ~ 12.8 x 10^*^ N/m^. The area of the unit cell in graphene is 



5.24 X 10 cm^ (a = 2.46 A). We then define the interlayer force constant 



per unit cell a = av2B with value 670 dyn/cm. In Ref. 



20l the rigid-layer shear frequency 



of graphite is obtained as the limit for ?sf — > oo of the GML frequency uj-^ and reads 
Uoo = 2A/«//i. Here /z = 2M/v2d = 7.6 x 10~® g/cm^ is the mass per unit area (M = 12 



u for C). Writing then Uc 



2a/ M, we compare with the expression of the rigid- layer 



shear mode in graphite co{E2g-^), Eq. (31) of Ref. 
expression in the form uj{E2g-^) = 



22 



It is straightforward to recast this 



2hxx/M, where h^x is a sum of interatomic vdW force 



221 we obtain 



constants between next-neighbor graphene planes. With the values of Ref. 
hxx = 654 dyn/cm, comparable with 670 dyn/cm 20|. For the graphene bilayer we obtain 
uJx^{J^ = 2) = \Jhxx/M and hence uj{E2g^) = y^iOXsC^ = 2). We recall that the values of 
these interlayer shear force constants have to be chosen ad hoc, they can not be obtained 



from currently accepted Lennard- Jones potentials for C-C vdW interactions [35|. Similar 
relations hold for the rigid-layer compression modes. We obtain uo{B2g^) = \J 2hzz/M and 
for the bilayer oJxX"^ = 2) = a;(i?2gi)v^- The force constant is given by h^z = 5755 dyn/cm. 



B. Relation with bulk dispersions 



We now establish quantitative relations between the shear-mode eigenfrequencies 
{u;a5(^)} of the multilayers GML and BNML and the TA and TO dispersions with wave vec- 
tor q = (0, 0, Qz) along A-F in graphite and h-BN, respectively. With the ?\f-layer system we 
associate a quantized wave length z/A^ = Kc/2 or equivalently a wave vector q'^(z^) = Atcu/'Nc, 
where u is an integer in the interval 1 < ly < 3\f/2. The center of mass of the multilayer 
system stays at rest in displacements associated with the wave vector We calculate 

the sum of phase factors of the corresponding displacement pattern. 
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which is zero since v is an integer in the interval [1, ?Nr/2]. Since the nearest distance between 
equivalent planes in the bulk materials is 2(c/2), we have the correspondence 

g.(-) ^ ^ = (15) 

We find that the !N — 1 rigid- layer shear modes {u;a5(?^)} of the Ji layer-system are obtained 
as the intersections of the two lowest phonon branches TA and TO in Fig. HJ^a) (red, full 
lines) with vertical lines located at g = (O, 0, ^^(z/)) , v integer G [1,^]- For even, the 
value V = \ corresponds to g = (0,0, ^), i.e. the A point of the 3D Brillouin zone. Hence 
one obtains the series of central points with constant frequencies in the fan diagrams. These 
are located at w = 30.4 cm~^ for the shear modes in GML [see Fig. EJ^a)], and at w = 38.5 
cm~^ for the shear modes in BNML [Fig. Et^b)]. The cuts of the graphite A-F branches 
resulting in the GML fan diagrams for K < 10 are shown in Fig. Hl^a) by vertical lines at 
given by Eq. f lTSj) : the values of = ^ as well as the corresponding X- values are labeled 
on the top horizontal axis. Note that the frequencies present for Ji are also present for all 
multiples of Ji. Also note that in the limit — > oo a cut of the graphite A-F branches 
at qz = 0, i.e. the F point, is reached, consistent with the observation that the frequency of 
the £'231 mode at the F point {u = 43.0 cm~^) corresponds to the limit value of the highest 
frequencies {wjj'^(!N)} of the lower fan diagram of Fig. (red, filled circles) mentioned 
before. On the other hand, the lowest frequencies {^^^(K)} of the lower fan diagram of 
Fig. [3]^a) (red, filled circles) evolve to zero for 'N — > 00, in agreement with the limit of the 
TA branch with g^ — y 0. We next consider the rigid-plane compression modes {a;Ac(^)} 
in GML, upper fan diagram of Fig. |3](a) (blue, open circles). Here too, the relation [Eq. 
f|T5|) ] with the LA and LO modes in graphite can be established. The two branches meet 
at the A point at 90.1 cm~^ which agrees with the value of the bilayer and rigid-plane 
compression eigenfrequency at gl = 0. The frequency a;(i?2gi) = 127.5 cm~^ of the LO 
branch at F in graphite corresponds to the limit for large !N of the sequence of highest 
eigenfrequencies of the compression modes in GML {u^^iJ^)} = {90.1 cm^^; 110.4 cm~^; 
117.8 cm~^; . . .; 125.9 cm~^; . . .} for K = {2; 3; 4; . . . ; 10; . . .}, respectively. Here the lowest 
frequencies {w^^([N")} tend to zero for K — t- 00, in agreement with the TA branch. The 
evolution of this mode with increasing number of layers has already been studied in one of 
our previous papers [2^. Since the -8231 mode in graphite is optically silent, this prediction 
has not been checked by experiment. However, from a group-theoretical analysis it has been 
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concluded [36| that in GML with Jsf even there are ?sf/2 compression modes with symmetry 
Aig that are Raman active. Also for the case of 'N uneven, infrared active modes should 
occur. Again, the eigenfrequencies plotted in the upper fan diagram in Fig. |3t^a) coincide 
with the intersections of the vertical lines q = (0, 0, Qz) with given by Eq. ( |T5l) with the 
LO and LA phonon branches along A-F in graphite. From a comparison of low-frequency 
out-of-plane phonon dispersions of [Nf-layer graphene at F and the low-frequency dispersions 



LO, LA of graphite along F-A, it has been inferred [37| that a relation like Eq. f|T5l) holds, 
however with u = 0,l,...,?\f — 1. This range of u overestimates the number of rigid plane 
compression modes by more than a factor 2. 

We now turn to BNML. Comparing with the low-frequency dispersions TA, TO and LA, 
LO of h-BN along A-F in the Brillouin zone, we find that the intersections obtained by 
means of Eq. f|T5|) determine again the eigenfrequencies of the multilayers. This holds as 
well for the shear modes {u;a5(?^)} as for the compression modes {w^cl^)} (lower and upper 
fan diagrams in Fig. [3]^b), respectively) and is illustrated in Fig. ID^b). 

Having established the relation between the A-F phonon branches TO, TA and LO, 
LA of the 3D materials and the eigenfrequencies {u;As(3\f)} and {u;a^(^)} at the F point 
of the multilayers for rigid-shear and rigid-compression modes [Eq. (|T5|l ]. it is possible to 
deduce u^N) master curves connecting the fan diagram frequencies (Fig. [3]). The following 
considerations hold for the shear modes as well as for the compression modes. In the former 
case ua and ur stand for uj^ and respectively, in the latter case for and Wp^. First, 
one needs to consider each fan diagram as a set of pairs of c<;(?sf) curves; at (K, co) = (2n, ua), 
with = 1, 2, 3, . . ., a pair of curves (one increasing, the other decreasing) originates. From 
Eq. (fT5|) it follows that the frequencies lying on the curves originating at {coa, 2n) are obtained 
by cutting the corresponding A-F branches of the 3D material at Qn = ^qz{n) = with 
!N = 2,3,4,.... The curves containing the highest and lowest frequencies (n = 1) are 
e.g. obtained by cuts at Qi = |, |, |, . . .. Secondly, the A-F 3D phonon branches can be 
extremely well approximated by second-degree curves. For the optical branches, satisfying 



dm 
dQ 



= 0, we put 

Q=0 



'^{Q = ^Qz = 0) = LOT, uj{Q = 1) = wa and 

"•^ y=u 

co+iQ) =ur-{uJr-ujA)Q^. (16) 
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For the acoustic branches, for which u^Q = 0) =0 and ijj{Q = 1) = coa, we assume 



-{Q)=ujAQ[{l + r)-rQ]^ 



(17) 



with r a dimensionless fit parameter. The superscripts '+' and '— ' refer to the increasing 
and the decreasing curves, respectively. Now inserting 



^ resuhs in 



(^) = - (wr - wa) 



(^) = 



2n' 

^1 + ^) 



(18) 
(19) 



The UJ^(y^^i) and w^(3^[) curves plotted in Figs.[3l^a) and[3](b) have been obtained by evaluating 
Eqs. (|T8l) and (|T9l) . For the latter, values of r = 0.165 and r = 0.166 were fitted for GML and 
BNML, respectively. The agreement between the multilayer phonon frequencies, obtained 
by diagonalizing the 6'N x 6K dynamical matrix, and the master curves, obtained by making 
cuts of the A-F 3D phonon branches, is perfect. Quadratic- form assumptions for the A-F 
rigid-shear and rigid-compression phonon branches obviously work extremely well (for both 
graphite and h-BN). Expression (|T71) for uj~{Q) is readily used to calculate the sound velocity 
of LA and TA phonons in the bulk materials. Substituting Q = ^q^ in Eq. ( 1T9|) we obtain 



V 



lim 



WA-(l + r). 

77 



(20) 



With c = 6.7 A, = 90.1 cm~^ and uijf^ = 30.4 cm~^ we obtain in case of graphite 
the longitudinal and transverse sound velocities Pla = 4.22 km/s and Vta = 1-42 km/s. 
The experimental values are 4.14(4) km/s and 1.48(6) km/s, respectively 38|. For h-BN, 
with c = 6.66 A, 00^^ = 86.3 cm^^ and ujjf^ = 38.6 cm~^, we obtain Vla = 4.02 km/s and 



Vta = 1-79 km/s, to be compared with the experimental values 39| 3.44(3) km/s and 1.84(6) 
km/s, respectively. Note further that the evolution of the rigid modes' highest frequencies 
(shear and compression) are given by 



cj+(K) = cor 



(wp-wa)^. 



(21) 



The procedure of obtaining the phonon eigenfrequencies of the multilayer system by 
making intersections of the A-F phonon branches of the 3D material is the analog of the 
zone-folding scheme where the phonon dispersion relations of ID carbon nanotubes are 
obtained from 2D graphene 40|, |41 |. 
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V. DISPLACEMENT CORRELATIONS 

Having determined the eigenfrequencies and eigenmodes for GML and BNML systems, 
we will calculate the temperature-dependent dynamic and static correlation functions. The 
knowledge of these functions is relevant for the interpretation of scattering experiments and 



is likely to be useful for friction experiments 13j . We use a quantum- mechanical formulation 



of lattice dynamics in the harmonic approximation and extend the standard theory of 3D 
crystals [27] to the case of multilayers. For a given ?\f-layer system we consider the time- 
dependent displacement operator Ui{n±, I, k; t). Here i G {x, y, z}, n± refers to the prismatic 
unit cell j2l|, / to the multilayer plane, < / < K — 1, k to the particle of mass m^, t stands 
for time. The expansion in terms of normal coordinates Qx{q±; t) reads 

Here X(nx, /, n) is the equilibrium position of the particle (n^, /, k) in the multilayer crystal, 
N± the number of unit cells. In terms of phonon creation (6^) and annihilation (6) operators 
the time- dependent normal coordinate reads 



One has the usual commutation relations for Bose operators 

[b,{q^),b{,i(f^)]=6.^y^6,y, (24) 
[bx{q±),byiq^)] = [bi{u),b[,{q'^)] = 0. (25) 

The thermal occupation of phonons with polarization A and frequency uJx{q±) is given by 

mq^)bxm - nx{q\) = , (26) 

where /3 = ksT, T is the temperature and k-Q the Boltzmann constant. 

In the case of rigid-layer displacements we retain at q± = those eigenmodes that satisfy 
Eqs. ([7]) - ffTOj) and we denote the eigenfrequencies ujxa{q± = 0) by ux^ and the occupation 
number n^^ (gl = 0) by nx^ . Here again Aq, = A^ refers to rigid-layer shear modes and Aq, = 
Ac to rigid-layer compression modes. The rigid-layer displacement-displacement dynamical 
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correlation function reads 

(u,(nx, /, t)ui{n^,l, 0))^ = M'^\Xa, 0)(gA.(0; t)Qx^{0; 0))^. 

(27) 

We have taken into account that only terms diagonal in Aq, contribute to the thermal average 
( ). Evaluation of the thermal average gives 

(QaJO; t)Q,^{0; 0))^ = K^^"^'^* + (1 + ^Aje"^-^-^*] . (28) 

Notice that the right-hand side of Eq. f l27|) is independent of n± (rigid layers!). In order to 
obtain the shear and compression correlation functions of the !N-layer system, we multiply 
both members of Eq. ( 1271) by the number of unit cells N±, sum over /, n and i G {x,y} or 
i = z. We call the result (ua(t)Ma(O))^ where a = s in case i G {x,y} and a = c for i = z. 
The result reads 

K(t)«a(0)>^ = /E(^^'^W^^'^(0)>3.' (29) 

Aq, 

where 

^ y 0)ef'^^(A., 0) d'^^A,, 0)d''^^(Ac, 0) 2 

f-^ 4^ m' 

with m = ^ {I5 2}. The result fl30l) is a consequence of Eqs. ([7]), ([8]) and ([9]). Note 

that the result is independent of the number of layers; one has / = 0.0806 u^^ or 0.0833 
u~^ for BNML or GNML, respectively (with m in atomic mass units u). The sum over 
on the right-hand side of Eq. (129|) depends on K. 

The static correlation functions are obtained by taking t = 0. From Eq. ( l28|l we get, by 
using Eq. (l26l) . 

,2\ fi /3^Ac 



and hence 

h f3hijj 
2~x 



-^(0)>. = /E^-^h^. (32) 



Ad 

Making use of relations (1181) and (1T9|) . we can write 



(33) 
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for 'N even {n = 1, 2, 3, . . . , j — 1) and 



«(0)>, = 9f E coth + coth j (34) 

for [N" odd (n = 1, 2, 3, . . . , ^^)- Here g = 2 accounts for the degeneracy in case of the shear 
modes while g = 1 in case of compression modes. We recall that Eqs. ( !33l) and ( !34l) apply 
to shear modes {A^} for ua and ur entering given by wj*-* and and to compression 
modes {Ac} for and Wp^. We have calculated the temperature-dependent mean-square 



displacements y (0))^^ for rigid-plane shear (a = s) and compression (a = c) modes by 
means of Eqs. ( 133|) and for GML and BNML. Results for a series of ?\f-layer systems are 
shown in Fig. [5l We recall that these results are obtained for rigid layers [Eq. Q], while the 
center of mass of the 3\f-layer system stays at rest. Hence the static displacement correlation 
functions (^^(O))^^ and (^^(O))^^ are a measure of the total amount of the relative shear and 
compression motion, respectively, between rigid layers. 

We next study the static displacement correlation function of the surface layer with label 
/ = 0. Instead of Eq. fl30|) . we have to consider 

/r (A.) = E (35) 

IK. 

where a = s for i G {x, y} and a = c for i = z. The static correlation function now reads 

(K°'(0))n = E/r(A<.)^coth^. (36) 

Act 

In Fig. [6] we show numerical results of the mean-square thermal displacements. 

From Figs. [5] and |6] we conclude that the average rigid-layer shear and compression dis- 
placements increase with increasing temperature and with layer number 'N. In Ref. [l3| the 
results of friction force microscopy experiments demonstrate that friction decreases monoton- 
ically with the number of layers. The mechanical origin for the observed effect is attributed 
to the fact that the sliding AFM tip causes out-of-phase deformations (puckering) of the 
surface sheet. The increased tip-sheet contact area or (and) the additional work required to 
move the puckered region forward lead to increased friction. This effect is more pronounced 
for thinner samples which exhibit a lower bending stiffness. On the other hand for thicker 
sheets the puckering is less prominent owing to the larger bending stiffness of the sheet [3]. 
Within this scenario it is suggested that some relative sliding between the topmost layer 
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FIG. 5: Temperature-dependent mean-square displacements y (n^(0))j^ for rigid-plane shear (a = 
s, top) and compression (a = c, bottom) for GML. The number of layers !N ranges from 2 (lowest 
curve) to 10 (upper curve). The results for BNML are very similar. 

and the material below occurs. This feature should increase the spacing of the stick-slip 
events. Our results (Fig. [6]) on the increase of the mean-square shear displacements of the 
surface layer with increasing !K are then compatible with the experimental findings that the 



spacing of the stick-slip events increases with increasing !N [13|. Concerning the increase of 
the vertical (compression) mean-square displacements with increasing K we are led to argue 
that those processes decrease the contact area between AFM tip and multilayer system and 
hence contribute to a decrease of friction with increasing [N. 

We close with a comment on dynamics. The Fourier transform of the time-dependent 
correlation function (^Ua{t)ua{0)') , a = s (c) is relevant for the interpretation of dynamic 
scattering laws. We define 



1 

2^ 



(37) 
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FIG. 6: Temperature-dependent mean-square displacements 



(0) 



(0))' 



for rigid-plane shear 



(a = s, top) and compression (a = c, bottom) for GML. The number of layers !N ranges from 2 
(lowest curve) to 10 (upper curve). The results for BNML are very similar. 



and obtain by means of Eqs. fl29l) and fl28|) 



{nxj{yj + waJ + (1 + ^aJ5(w - waJ] • 



(38) 



Here huj stands for the energy transfer of the scattering particle (photon or neutron) to the !N- 
layer system. The first term within the square bracket represents an energy absorption by the 
scattering particle (anti-Stokes process) and the second term an energy loss (Stokes process) 
which becomes dominant at low T. Expression f l38|) comprises all shear or compression 
motion resonances of a given Jsf-layer system. So far the highest value shear resonances {w^^^ } 
have been detected by experiment 2C|. While these experiments have been carried out at 
room temperature, it might be necessary to go to lower T in order to detect the resonances 
at lower frequencies. Also the Raman-active compression modes in even !N multilayers 
symmetry Aig, are a challenge for further experiments. 



22 



VI. CONCLUDING REMARKS 

We have given a theoretical investigation of the low-frequency phonon dispersions in 
crystalline layered materials. These phonons, associated with rigid-plane motions, show 
universal behavior which applies to metallic (GML) as well as to ionic, insulating (BNML) 
systems. The frequency spectra have been represented in the form of fan diagrams for 
compression (also called stretching) and shearing motions. For a system of "N layers one 
distinguishes !N — 1 compression modes and N — 1 doubly degenerate shear modes with 
frequencies {(X'Ac(?^)} and {a;A^(?^)}, respectively. The fan diagrams (see Fig. [3]) are centered 
around a series of frequency points given by the bilayer frequencies oJx^l'N = 2) and ujx^CN = 
2) appearing for systems with an even number of layers 'N. The fan diagram associated with 
compression is centered around higher frequencies than the fan diagram associated with shear 
motion in both GML and BNML. For both shearing and compression the sequences of highest 
frequencies {a;^_^(IN")} and {u'^^{'N)} have as limits for K — > oo the bulk material frequencies 
uj{E2gJ and uj{B2g^) respectively at the F point of the 3D Brillouin zone. In case of GML the 
series of shear frequencies {w^ (^)} up to K = 11 has been measured by Raman scattering 
[2^. On the other hand the sequences of lowest frequencies {a;^^_^([N")} and {a;5^^(3\f)} have 
limit values for K — > 00. Comparison with the low-frequency dispersions along the F- 
A line in the Brillouin zone of hexagonal layered 3D materials shows that the frequencies 
uxX^ = 2) and ujx^^ = 2) of the bilayer agree with the bulk frequencies uiTO) and w(LO), 
respectively, at the A point. In addition one has the relations ^/2ux^{J^ = 2) = uj{B2g^) and 
\/2ujxX^ = 2) = uj{E2gj). These relations are a consequence of the fact that the interlayer 
force constants hzz for compression and hxx for shearing are only effective between next- 
neighbor rigid planes. Interactions between more distant rigid planes are negligible. Note 
rigid-layer shear modes in GML this conclusion was drawn from Raman 



that for the case o: 

scattering results [20|. We have attributed the absence of longer distance interactions to 
screening effects. In GML the screening is due to the metal nature (vr electrons), in BNML 
the screening is due to the overall charge neutrality and the plane rigidity. 

We further have explored the relations between the LO and LA phonons of the bulk 
materials along q = [0, 0, g^] and the fan diagram frequencies for compression modes in 
!N-layer systems; similar relations exist between TO and TA phonons along A-F and the 
fan diagram frequencies for shearing modes. In both cases the 3\f — 1 rigid-layer frequencies 
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{wac(^)} and {cjas(^)} ^ire obtained as intersections of the LO, TA and TO, TA phonon 
branches, respectively, with vertical lines at discrete positions ^^(z/) along F-A. We have 
obtained master curves which allow to derive the fan diagrams of GML and BNML for any 
given [N". 

Finally we have calculated static and dynamic correlation functions for rigid-plane mo- 
tions. We have studied correlations as functions of T and of !N. Our results, which exhibit 
again a large similarity between GML and BNML, might be of relevance for the understanc- 



ing on an atomistic level of the results of force friction experiments on thin-layer sheets 13 |. 
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